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ABSTRACT 

A description is given for preserving V • B = in a magnetohydrodynamic (MHD) 
code that employs the upwind, Total Variation Diminishing (TVD) scheme and the 
Strang-type operator splitting for multi-dimensionality. The method is based on the 
staggered mesh technique to constrain the transport of magnetic field: the magnetic 
field components are defined at grid interfaces with their advective fluxes on grid edges, 
while other quantities are defined at grid centers. The magnetic field at grid centers 
for the upwind step is calculated by interpolating the values from grid interfaces. The 
advective fluxes on grid edges for the magnetic field evolution are calculated from the 
upwind fluxes at grid interfaces. Then, the magnetic field can be maintained with 
V • B = exactly, if this is so initially, while the upwind scheme is used for the update of 
fluid quantities. The correctness of the code is demonstrated through tests comparing 
numerical solutions either with analytic solutions or with numerical solutions from 
the code using an explicit divergence-cleaning method. Also the robustness is shown 
through tests involving realistic astrophysical problems. 

Subject headings: magnetohydrodynamics: MHD - methods: numerical 
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1. INTRODUCTION 

The current bloom in computational astrophysics has been fed not only by dramatic 
advances in computer hardware, but also by comparable developments in improved algorithms. 
Nowhere has that been more important than in methods to solve the equations of compressible 
magnetohydrodynamics (MHD). That system is the most applicable to describe a vast array 
of central astrophysical problems. But, MHD has presented a special challenge, because of the 
complexity presented by three non-isotropically propagating wave families with wide ranging 
relative characteristic speeds, and the associated need to solve a reduced set of Maxwell's equations 
along with the equations of compressible continuum fluid dynamics. 

The condition V • B = is a necessary initial constraint in multi-dimensional MHD flows 
and should be preserved during their evolution. While the differential magnetic induction 
equation formally assumes V • B = 0, non-zero V • B can be induced over time in numerical 
simulations by numerical errors due to discretization and operator splitting. This is because, 
even though conventional numerical schemes may be exactly conservative of the advective fluxes 
in the induction equation, nothing maintains the magnetic fluxes in the sense of Gauss's Law. 
That is, nothing forces conservation of zero magnetic charge within a finite cell during a time 
step. Numerical non-zero V • B usually grows exponentially, causing an anomalous force parallel 
to the magnetic field and destroying the correct dynamics of flows, as pointed by Brackbill & 
Barnes (1980). Those authors, along with Zachary et al. (1994) show that the use of a modified, 
non-conservative form of the momentum equation can keep the non-zero V • B small enough that 
no further correction is necessary for this purpose. However, the modified form is not suitable for 
some schemes, and more importantly may result in unphysical results due the non-conservation of 
momentum (see, e.g., LeVeque 1997). 

Several methods have been suggested and used to maintain V • B = in MHD codes. We 
mention four here. In the first method, vector potential is used instead of magnetic field in the 
induction equation (see e.g., Clarke et al. 1986; Lind et al. 1989). Although V • B = is ensured 
through the combination of divergence and curl operations, the method results in the second-order 
derivatives of the vector potential in the Lorentz force term of the momentum equation. So, in 
order to keep second-order accuracy, for instance, the use of a third-order scheme for spatial 
derivatives is required (for detailed discussion, see Evans h Hawley 1988 and references therein). 
In the second method, the MHD equations are modified by adding source terms and any non-zero 
V ■ B is advected away from the dynamical region (see Powell 1994 for details). That method 
works well for some problems with open boundaries but not in others, including those with 
periodic boundaries. In the third method, an explicit divergence-cleaning scheme is added as a 
correction after the step to update fluid quantities (see e.g., Zachary et al. 1994; Ryu et al. 1995a; 
Ryu et al. 1995b). The method works well if boundary effects are negligible or the computational 
domain is periodic and if the grid used is more or less regular. Otherwise, however, the scheme is 
not easily adaptable. In the fourth method, the transport of magnetic field is constrained by the 
use of a staggered mesh: some quantities including magnetic field components are defined on grid 
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interfaces while other quantities are defined at grid centers. The method has been successfully 
implemented in schemes based on an artificial viscosity (see e.g., Evans & Hawley 1988; DeVore 
1991; Stone & Norman 1992). However, since it is "unnatural" to stagger fluid quantities in 
Riemann-solver-based schemes, that approach has only recently been applied successfully in such 
schemes. 

Conservative, Riemann-solver-based schemes, which are inherently upwind, have proven to 
be very effective for solving MHD equations as well as hydrodynamic equations. These schemes 
conservatively update the zone-averaged or grid-centered fluid and magnetic field states based on 
estimated advective fluxes of mass, momentum, energy and magnetic field at grid interfaces using 
solutions to the Riemann problem at each interface. MHD examples include Brio & Wu (1988), 
Zachary & Colella (1992), Zachary, et al. (1994), Dai & Woodward (1994a, 1994b), Powell (1994), 
Ryu & Jones (1995), Ryu, et al. (1995a), Powell et al. (1995), Roe & Balsara (1996), Balsara 
(1997), and Kim et al. (1998). Brio & Wu applied Roe's approach to the MHD equations. Zachary 
and collaborators used the BCT scheme to estimate fluxes, while Dai & Woodward applied the 
PPM scheme to MHD. Ryu and collaborators extended Harten's TVD scheme to MHD. Powell 
and collaborators developed a Roe-type Riemann solver with an eight-wave structure for MHD 
(one more than the usual number of characteristic MHD waves), one of which is used to remove 
non-zero V • B. Balsara used also the TVD scheme to build an MHD code. 

The upwind schemes share an ability to sharply and cleanly define fluid discontinuities, 
especially shocks, and exhibit a robustness that makes them broadly applicable. But, due to 
the feature of the upwind schemes that zone-averaged or grid-centered quantities are used to 
estimate fluxes at grid interfaces, the staggered mesh technique has been slow to be incorporated 
for magnetic flux conservation. Instead, the explicit divergence-cleaning scheme has been more 
commonly used. However, recently Dai & Woodward (1998) suggested an approach to incorporate 
the staggered mesh technique in the upwind schemes. It relies on the separation of the update 
of magnetic field from that of other quantities. Quantities other than magnetic field are updated 
in the upwind step either by a split or an un-split method. Then the magnetic field update is 
done through an un-split operation after the upwind dynamical step. Magnetic field components 
are defined on grid interfaces, and their advective fluxes are calculated on grid edges using the 
time-averaged magnetic field components at grid interfaces and the other time-averaged quantities 
at grid centers through a simple spatial averaging. For the upwind dynamical step, the values of 
magnetic field at grid centers are interpolated from those on grid interfaces. 

In this paper, we describe an implementation of the Dai & Woodward approach into a 
previously published upwind MHD code based on the Harten's Total Variation Diminishing (TVD) 
scheme (Harten 1983) which the present authors developed (Ryu & Jones 1995; Ryu, et al. 1995a; 
Ryu, et al. 1995b; Kim et al. 1998). The TVD scheme is a second-order-accurate extension of the 
Roe-type upwind scheme. The previous code employed an explicit divergence-cleaning technique in 
multi-dimensional versions, and has been applied to a variety of astrophysical problems including 
the MHD Kelvin-Helmholtz instability (Frank et al. 1996; Jones et al. 1997), the propagation of 
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supersonic clouds (Jones et al. 1996), and MHD jets (Frank et al. 1998). However, the range of 
application for the code has been limited due to the restrictions on boundary condition and grid 
structures, as pointed out above. In this paper we address those limitations by incorporating the 
staggered mesh algorithm to keep V ■ B = 0. However, instead of calculating the advective fluxes 
for the magnetic field update as Dai & Woodward (1998), we calculate the fluxes at grid edges 
using the fluxes at grid interfaces from the upwind step. The advantage of our implementation is 
that the calculated advective fluxes keep the upwindness in a more obvious way. We show that 
our new code performs at least as well as the previous version in direct comparisons, but also that 
it effectively handles problems that could not be addressed with the original code. We intend this 
paper to serve as a reference for works which use the code for astrophysical applications. In §2, the 
numerical method is described, while several tests are presented in §3. A brief discussion follows 
in 84. 



2. IMPLEMENTATION OF THE DIVERGENCE-FREE STEP 

We describe the magnetic field update step in two-dimensional plane-parallel geometry. 
Extensions to three-dimension and other geometries are trivial. The induction equation in the 
limit of negligible electrical resistivity is written in conservative form as 

dB d 

~df + dy' {BxVy ~ ByVx) = °' (1) 

and 

dB d 

-g L + ^(B y v x -B x v y ) = 0. (2) 

The full MHD equations in conservative form can be found, for instance, in Ryu et al. (1995a). In 
typical upwind schemes, including the TVD scheme for MHD equations as well as hydrodynamic 
equations, fluid quantities are zone-averages or defined at grid centers. Their advective fluxes are 
calculated on grid interfaces through approximate solutions to the Riemann problem there, based 
on interpolated values. 

Here, we define the magnetic field components on grid interfaces, b x i j and b yi ij, while all 
the other fluid quantities are still defined at grid centers (See Fig. 1 for the notations used in this 
paper). For use in the step of calculating the advective fluxes by the TVD scheme, the magnetic 
field components at grid centers, which are intermediate variables, are interpolated as 

Bx,i,j = 2 few + b x ,i-i,j) > (3) 

and 

By,i,3 = 2 $y,i,j + b y ,i,j-i) ■ (4) 

Since the MHD code based on the TVD scheme has second-order accuracy, the above second-order 
interpolation should be adequate. If non-uniform grids are used, an appropriate interpolation 
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of second-order should be used. With the fluid quantities, including these magnetic field values, 
given at grid centers, the TVD advective fluxes are used to update the fluid quantities from the 
time step n to n + 1 

(5) 



<f +1 



T T n* 1 

±J y 1 -'xHi ) j-i 



as described in Ryu & Jones (1995) and Ryu et al. (1995a). Strang- type operator splitting is used, 
so that the operation L x L y is applied from n + 1 to n + 2. Here, q is the state vector of fluid 
variables. 

The advective fluxes used to update the magnetic field components at grid interfaces are 
calculated also during the TVD step from the MHD Riemann solution with little additional cost. 
In the x-path the following flux is computed 



- 1 / \ Aa; 7 

fx,i,j = g \ B y,i,j V x,i,j + B y,i+l,j v x,i+l,j) ~ 2~A^T $k,i+\ ,j R k,i+± j( 5 )' 

k=l 



(6) 



and in the y-path the following flux is computed 



fv>h3 - 2 { B x,i,j v y,i,j + B x,i,j+l v y,i,j+l) 2 At n ^ 

k=l 



(7) 



Here, 3 h ■ , i ■ and j3, ■ ■ , i are the quantities computed at grid interfaces. As described in Ryu & 
Jones (1995), (3 ki+ ij used in the the x-path is calculated as follows: 
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a k,i+lj ~ ^kj+lj • Wi+lJ 



0, for ciju „• , i „• = 



g k>iJ = sign^ ^i^.) max [o, min { 
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Ax kji+^j Ax kji+^j 



a kM\,3' 



(8) 
(9) 

(10) 

(11) 
(12) 

(13) 



|X|, for |x| > 2e fc _ 

Pkii+± use< ^ i n the y-path is calculated similarly. R? i .(5) and i?" . i (5) are the fifth 

components ("B y " and U B X ", respectively for the two passes) of the right eigenvectors and L^'s 
are the left eigenvectors. They are computed on grid interfaces and given in Ryu & Jones (1995). 
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a£'s are the speeds of seven characteristic waves which are also computed on grid interfaces. Along 
the x-path, they are in non-increasing order 

ai,7 = v x ±c f , a 2fi = v x ±c a , a 3j5 = v x ±c s , a 4 = v x . (14) 

a^'s along the y-path is computed by replacing v x with v y . Here, c/, c a , and c s are local fast, 
Alfven, and slow speeds, respectively, e^s are the internal parameters to control dissipation in 
each characteristic wave and should be between and 0.5 (see the next section). The time step 
At n is restricted by the usual Courant condition for the stability. 

Using the above fluxes at grid interfaces, the advective fluxes, or the z-component of the 
electric field, on grid edges (see Fig. 1 for definition) are calculated by a simple arithmetic average, 
which still keeps second-order accuracy: namely, 



Then, the magnetic field components are updated as 

At n 

Ay 

and 



b ti = K,, - — (flij - (kj ;) , (16) 



b vS = K,, + 37 («« - fii-u) • (17) 

Note that the VL terms include information from seven characteristic waves. It is also clear that 
the net magnetic flux across grid interfaces is exactly kept to be zero at the step n + 1 

j s b n+1 ■ dS = (6£+} - b^_ ld )Ay + (6J+J - b^Ax = 0, (18) 

if it is zero at the step n. 

The reason why we take the fluxes in (6) and (7) from the upwind fluxes for the transport of 
the magnetic field at grid centers is this. As can be seen in equation (2) along the x-path, it is 
d(B y v x )/dx that contains the advective term and requires modification of fluxes to avoid numerical 
problems. d(B x v y )/dx causes no problems. The same argument is applied to d(B x v y ) / dy along 
the y-path. We note that in the above scheme, the results of one-dimensional problems calculated 
with the two-dimensional code reduce to those calculated with the one-dimensional code, as it 
should be. For instance, in shock tube problems with structures propagating along a coordinate 
axis, the two-dimensional code produces exactly the same results as those given in Ryu & Jones 
(1995). However, with Dai & Woodward's (1998) advective fluxes, that is not necessarily true. 

The code runs at ~ 400 Mflops on a Cray C90, similar to the previous code (Ryu et 
al. 1995a). It corresponds to an update rate of ~ 1.2 x 10 5 zones s _1 for the two-dimensional 
version, about a 20% speed-up compared to the previous code, due to the absence of the explicit 
divergence-cleaning step. 
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3. NUMERICAL TESTS 

The numerical scheme described in the last section was tested with two-dimensional problems 
in plane-parallel and cylindrical geometries in order to demonstrate its correctness and accuracy 
as well as to show its robustness and flexibility. In all the tests shown, we used the adiabatic 
index 7 = 5/3 and a Courant constant C cour = 0.8. For the internal parameters, £fc's, of the TVD 
scheme (Ryu & Jones 1995; Ryu et al. 1995a), £17 = 0.1 — 0.2 (for fast mode), £2,6 = 0.05 — 0.1 
(for slow mode), £3,5 = — 0.05 (for Alfven mode), and £4 = — 0.1 (for entropy mode) were used. 
However, the test results are mostly not very sensitive to C cour and £& values. 

3.1. Shock Tube Problems 

We first tested the code with MHD shock tube problems placed diagonally on a two- 
dimensional plane-parallel grid. The correctness and accuracy are demonstrated through the 
comparison of the numerical solutions with the exact analytic solutions from the non-linear 
Riemann solver described in Ryu & Jones (1995). The calculations were done in a box of x = [0, 1] 
and y = [0, 1], where structures propagate along the diagonal line joining (0,0) and (1, 1). Two 
examples are presented. The first, shown in Fig. 2a, includes only two (x and y) components of 
magnetic field and velocity, so that they are confined in the computational plane. The second, 
shown in Fig. 2b, includes all three vector field components. The numerical solutions are marked 
with dots, and the exact analytic solutions are drawn with lines. Structures are measured along 
the diagonal line joining (0,0) and (1, 1). The plotted quantities are density, gas pressure, total 
energy, v» (velocity parallel to the diagonal line; i.e., parallel to the wave normal), v± (velocity 
perpendicular to the diagonal line but still in the computational plane), v z (velocity in the 
direction out of plane), and the analogous magnetic field components, B», B±, and B z . 

In Fig, 2a, the initial left state is (p, v±, v z , B±, B z , E) = (1, 10, 0, 0, 5/y/4n, 0, 20) 
and the initial right state is (1, —10, 0, 0, 5/\/47r, 0, 1), with B\\ = 5/\/47r. The calculation 
was done using 256 x 256 cells, and plots correspond to time t = 0.08\/2. The structures are 
bounded by a left and right facing fast shock pair. There are also a left facing slow rarefaction, a 
right facing slow shock and a contact discontinuity. All are correctly reproduced. The captured 
shocks and contact discontinuity here are very similar to those with the code using an explicit 
divergence-cleaning scheme, which was shown in Fig. 1 of Ryu et al. (1995a). 

In Fig. 2b, the initial left state is (p, v\\, v±, v z , B±, B z , E) = (1.08, 1.2, 0.01, 0.5, 3.6/a/47t, 
2/a/47T, 0.95) and the initial right state is (1, 0, 0, 0, 4/v / 4vr, 2/ v / 4vr, 1), with B\\ = 2/ v / 4vr. 
Again the calculation was done using 256 x 256 cells, and plots correspond to time t = 0.2\/2. 
Fast shocks, rotational discontinuities, and slow shocks propagate from each side of the contact 
discontinuity, all of which are correctly reproduced. Again, the structures are captured very 
similarly to what we found with the code using an explicit divergence-cleaning scheme shown in 
Fig. 2 of Ryu et al. (1995a). 
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3.2. The Orszag-Tang Vortex 

As a second, and truly multi-dimensional test, we followed the formation of the compressible 
Orszag-Tang vortex. The problem was originally studied by Orszag & Tang (1979) in the context 
of incompressible MHD turbulence, and later used to test compressible MHD codes by Zachary, et 
al. (1994), Ryu et al. (1995a), and Dai & Woodward (1998). Comparison of our new solution with 
previous ones demonstrates the correctness of the new code in this problem. 

The test was set up in a two-dimensional periodic box of x = [0, 1] and y = [0, 1] with 
256 x 256 cells. Initially, velocity is given as v = v D [— sm(2iry)x + sin(2-7rx)y] and magnetic field 
as B = B Q [— sin(27T?/)f + sin(47rx)y] with v = 1 and B Q = Uniform background density 

and pressure were assumed with values fixed by M 2 = v / {iPo/ Po) = 1, /? = p /(B 2 /2) = 10/3 
and 7 = 5/3. 

Fig. 3 shows the gray scale images of gas pressure, magnetic pressure, compression, V • v, and 
vorticity, (V x v) z at time t = 0.48, as well as the line cuts of gas pressure and magnetic pressure 
through y = 0.4277. The structures, including fine details, match exactly with those in Ryu et 
al. (1995a), proving the correctness of our code. Although only approximately the same initial 
conditions and epoch as those in Zachary, Malagoli, & Colella (1994) and Dai & Woodward (1998) 
are used, the images show that the overall shape and dynamics match closely with those solutions, 
as well. 

3.3. Propagation of a Supersonic Cloud 

As an initial test for the robustness and flexibility of the code in a practical, astrophysical 
application, we simulated supersonic cloud propagation through a magnetized medium. We 
present three simulations differing in initial magnetic field orientation. The first two reproduce 
Models A2 and T2 from Jones et al. (1996), in order to compare with previous calculations. The 
third is a new simulation. All three were computed on a Cartesian grid. In the first two cases the 
magnetic field is successively parallel to the cloud motion (aligned case - A2) and perpendicular 
to it (transverse case - T2). In the third, we present a new case with the magnetic field making an 
angle 6 = 45° with the cloud velocity (i.e., an oblique case). 

For these calculations we used the same physical parameters as in Jones et al. (1996). The 
cloud is initially in pressure equilibrium with the background medium and p Q = I/7 throughout. 
In addition, p am bient = 1) and the cloud density p c = X Pambient = 10- A thin boundary layer with 
hyperbolic tangent density profile and characteristic width of two zones separates the cloud from 
the background gas. The background sound speed is a am Uent = {iPol ' P ambient) 1 ^ 2 = 1- At the 
outset of each simulation, the ambient gas is set into motion around the cloud with a Mach number 
M = 10. The magnetic field lies in the computational plane and is initially uniform throughout it. 
Its strength corresponds to (3q = 4. Therefore, it is (B x , B y , B z ) = (0.55,0,0) for the aligned case, 
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(B x , By, B z ) = (0,0.55,0) for the transverse case and (B x , B y , B z ) = (0.39,0.39,0) for the oblique 
case. As in Jones et al. (1996) for the aligned and the transverse field cases the computational 
domain is [x,y] = [10,5], whereas for the oblique field case it is [x,y] = [20, 10]. The resolution is 
always 50 zones per initial cloud radius, R c ioud = 1- 

Images of the aligned and transverse cases are presented in Fig. 4a. Left panels correspond 
to the transverse case, and show density images (top two panels) and magnetic field lines (bottom 
panels) for two evolutionary times, namely t/t& c = 2, 6, where t& c is the bullet crushing time (see 
Jones et al. 1996 for detail). These are approximately the same times as those shown in Figs. 1 
and 2 in Jones et al. (1996). Figures representing the aligned field case are analogously illustrated 
in the right panels. As we can see there is a general agreement in both the density distribution 
and the magnetic field structure between the cases in Fig. 4a and the corresponding cases (T2 
and A2) in Jones et al. (1996). Minor differences appear in the details of the cloud shape and 
the magnetic field adjacent to the cloud for the aligned field case at t = 6tf, c . As pointed out in 
Jones et al. (1996), the nonlinear evolution of these clouds depends very sensitively on the exact 
initial perturbations and their growth. For both sets of simulations the perturbations develop 
out of geometrical mismatches between the cloud and the grid. We showed in that paper, for 
example, that consequently a simple shift of the initial cloud center by 0.5 zone on the x axis 
caused differences much greater than those illustrated here. Similarly, even minor changes in the 
field adjacent to the cloud near the start of the calculation or in the dissipation constants used 
can be expected to lead to observable changes in the detailed cloud features. Thus, by considering 
different schemes to keep V • B = used as well as different values of C cour and e^'s used in 
the different sets of simulations, we judge the agreement between the two codes to be good. In 
addition, the new code seems better able to handle the extreme rarefaction that forms to the 
rear of the cloud immediately after it is set in motion. That is a severe test, since the plasma (3 
abruptly drops from values larger than unity to values smaller than (3 ~ 10~ 2 . 

The simulation of a cloud interacting with an oblique magnetic field offers a good example of 
the increased flexibility of the new code. This situation is more realistic than the other two, but 
difficult to simulate with the old code because of its lack of a suitable periodic space for solving 
Poisson's equation in the explicit divergence-cleaning step. Since the aligned and transverse 
field cases differ considerably in their dynamics, it is astrophysically important to be able to 
investigate the general case of an oblique magnetic field. Fig. 4b illustrates the properties of one 
such simulation with the new code. Further details are discussed in Miniati et al. (1998b). Top 
and bottom panels correspond to density distribution and field line geometry respectively for 
two different evolutionary times (again t/% c = 2, 6). As we can see the evolution of the oblique 
case produces several features analogous to the previous transverse field case. In particular the 
magnetic field lines drape around the cloud nose and form an intense magnetic region there, due 
to field line stretching. In this fashion the field lines compress the already shock crushed cloud 
and prevent the rapid growth of the Kelvin-Helmholtz and Rayleigh- Taylor instabilities (Jones et 
al. 1996). However, the broken symmetry across the motion axis also generates uneven magnetic 
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field tension that causes some rotation and lateral motion of the cloud. In addition, it enhances 
turbulent motions in the wake and, therefore, the onset of the tearing mode instability and 
magnetic reconnection there. 



3.4. Jets 

As a final test to confirm the robustness of the new code and to demonstrate its application 
with a different grid geometry, we illustrate the simulation of a light cylindrical MHD jet with a 
top-hat velocity profile. The jet enters a cylindrical box of r = [0, 1] and z = [0, 6.64] at z = 0. 
The grid of the box is uniform with 256 x 1700 cells and the jet has a radius, rj e t, of 30 cells. 
The ambient medium has sound speed a am bient = 1 arid poloidal magnetic field (B^ = B r = 0, 
B z = B am i,ient) with magnetic pressure 1% of gas pressure (plasma Pambient = 100). The jet 
has Mach number M jet = v jet /a ambient = 20, gas density contrast p jet / Pambient = 0.1, and gas 
pressure in equilibrium with that of the ambient medium. It carries a helical magnetic field with 
B T = 0, B^ = 2 x B ambient (r/r jet ), and B z = B ambie nt- The jet is slightly over-pressured due to 
the additional B^ component, but the additional pressure is too small to have any significant 
dynamical consequences. The simulation was stopped at t = 2.2 when the bow shock reached the 
right boundary. It takes about about 20 CPU hours on a Cray C90 processor or about 90 CPU 
hours on a SGI Octane with a 195 MHz R10000. 

Fig. 5 shows the images of the log of the gas density and total magnetic field pressure 
(Fig. 5a) and the velocity vectors and the r and z magnetic field vector components (Fig. 5b) 
at five different epochs, t = 0.3, 0.8, 1.3 1.8, and 2.2. The length of velocity arrows is scaled as 
vTu| and that of magnetic field arrows as 

IB] 1 / 4 . Thefi gures exhibit clearly the complexity and 
unsteadiness of the flows. By viewing an animation of the simulation, which is posted in the web 
site |ht t p : / /canopus . chungnam . ac . kr / ryu /test j et / test j et . html , it becomes obvious that all of the 



structures are ephemeral and/or highly variable. The most noticeable structures are the bow shock 
of the ambient medium and the terminal shock of the jet material. In addition, the jet material 
expands and then refocuses alternately as it flows and creates several internal oblique shocks, as 
described in many previous works, for example, Lind et al. (1989). The terminal and oblique 
shocks are neither steady nor stationary structures. The oblique shocks interact episodically with 
the terminal shock, resulting in disruption and reformation of the terminal shocks. The terminal 
shock includes a Mach stem, so that the jet material near the outside of the jet exits through 
the oblique portion of the shock. That material carries vorticity and forms a cocoon around the 
jet. The vorticity is further developed into complicated turbulent flows in the jet boundary layer, 
which is subject to the Kelvin-Helmholtz instability. There are distinct episodes of strong vortex 
shedding which coincide with disruption and reformation of the terminal shock. Its remnants are 
visible as rolls in the figures. 

Although the total magnetic field in the back- flow region is strong compared to B amb i ent) 
as can be seen in the Pf, images, the components B r and B z are comparatively small, as can be 
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seen in the vector plots. This is because reconnection induced by the complicated turbulent flow 
motion of the jet material has frequently annihilated B r and B z , at the same time that has 
been enhanced by stretching. In an axis-symmetric calculation, the B^ component cannot be 
modified by reconnection, since it is decoupled from the other two magnetic field components. We 
emphasize that the details of the magnetic field configuration are sensitive to the assumed helical 
field within the incoming jet, so that our test results are representative only. 

Good agreement of this simulation with previous works such as Lind et al. (1989) provides 
another confirmation of validity and applicability of the new code. Detailed discussion of 
comparable jet simulations carried out with this new code in the context of radio galaxies, 
including acceleration and transport of relativistic electrons, have been reported in Jones et 
al. (1998). 

4. DISCUSSION 

For ordinary gasdynamics development of conservative, high order, monotonicity-preserving, 
Riemann-solution based algorithms, such as the TVD scheme employed here, provided a key 
step by enabling stable, accurate and sharp capture of strong discontinuities expected in 
compressible flows, while efficiently following smooth flows with a good economy of grid cells. The 
methods maintain exact mass, energy and momentum conservation and seem to do a good job 
of representing sub-grid-scale dissipation processes {e.g., Porter and Woodward 1994). Recent 
extension of those methods to MHD have also shown great promise, since they offer the same 
principal advantages as in gasdynamics. The main disadvantage of the Riemann methods in MHD 
were, until now, that they are basically finite volume schemes, so that they depend on knowing 
information averaged over a zone volume, or equivalently in second order schemes, at grid centers. 
The problem this presented came from the fact that the conservation of magnetic charge depends 
on a surface integral constraint, which is not guaranteed by the conservation of advective fluxes 
used in the remaining set of MHD relations. As discussed in the introduction this can lead to 
physically spurious results. 

Consequently, it is a significant advance to develop accurate, efficient and robust schemes for 
maintaining zero magnetic charge that are adaptable to Riemann-based methods. The method 
discussed in this paper seems to be an excellent choice. Since it exactly conserves the surface 
integral of magnetic flux over a cell and does it in an upwind fashion, it represents a class of 
techniques that have come to be called "Method of Characteristics, Constrained Transport" 
or "MoCCT". In this paper we outline a specific implementation of this scheme inside a 
multi-dimensional MHD extension of Harten's TVD scheme. With our prescription it should be 
straightforward for other workers to accomplish the same outcome. Through a varied bank of test 
problems we have been able to demonstrate the accuracy and the flexibility of the methods we 
have employed. Thus, we believe this code and others like it offer great potential for exploration of 
a wide variety of important astrophysical problems. Already the code described in the paper has 
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been used successfully in Jones et al. (1998), Miniati et al. (1998a), and Miniati et al. (1998b) to 
study propagation of cylindrical MHD jets including the acceleration and transport of relativistic 
electrons, and to study the propagation and collision between interstellar plasma clouds. 
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Fig. 1. — Notations for the flow variables used in the paper. The centered magnetic field, B, and 
the velocity, v, are defined at grid centers. The faced magnetic field, b, and the upwind fluxes, /, 
are defined on grid interfaces. The advective fluxes for the magnetic field update, Q, is defined at 
grid edges. 

Fig. 2a. — Two-dimensional MHD shock tube test. Structures propagate diagonally along the line 
from (0,0) to (1, 1) in the x — y plane. The initial left state is (p, v\\, v±, v z , B±, B z , E) = (1, 10, 
0, 0, 5/a/4tF, 0, 20) and the initial right state is (1, -10, 0, 0, 5/\/4tt, 0, 1) with B\\ = 5/VAtt (same 
test as Fig. 1 in Ryu et al. 1995a). The plots are shown at time t = 0.08\/2 along the diagonal 
line joining (0,0) and (1, 1). Dots represent the numerical solution from the code described in §2 
with 256 x 256 cells. Lines represent the exact analytic solution from the nonlinear Riemann solver 
described in Ryu & Jones (1995). 

Fig. 2b. — Two-dimensional MHD shock tube test. Structures propagate diagonally along the line 
from (0,0) to (1, 1) in the x — y plane. The initial left state is (p, v\\, v±, v z , B±, B z , E) = (1.08, 
1.2, 0.01, 0.5, 3.6/a/4tF, 2/V47T, 0.95) and the initial right state is (1, 0, 0, 0, A/V&r, 2/v / 4vf, 1) 
with Bii = 2/\/4vr (same test as Fig. 2 in Ryu et al. 1995a). The plots are shown at time t = 0.2\/2 
along the diagonal line joining (0,0) to (1, 1). Dots represent the numerical solution from the code 
described in §2 with 256 x 256 cells. Lines represent the exact analytic solution from the nonlinear 
Riemann solver described in Ryu & Jones (1995). 

Fig. 3. — Gray scale images of gas pressure (upper left), magnetic pressure (upper right), V-v (lower 
left), and (V x v) z (lower right) in the compressible Orszag-Tang vortex test. White represents 
high (or positive) values and black represents low (or negative) values. The calculation has been 
done in a periodic box of x = [0, 1] and y = [0, 1] with 256 x 256 cells. The initial configuration 
is p = 25/367T, p = 5/127T, v = — sin(27ry)x + sin(27rx)y, and B = [— sin(27ry)x + sin(47rx)y] /V^tt, 
and the images shown are at t = 0.48. The line plots show the profiles of gas pressure and magnetic 
pressure along y = 0.4277. 

Fig. 4a. — A supersonic cloud moving through a magnetized medium and computed on a Cartesian 
grid. The initial cloud radius in each case is 50 cells. Upper panels show logarithmic gray scale 
images of gas density, with white referring to high density values and black to low ones. Lower 
panels illustrate magnetic field lines obtained as contours of the magnetic flux function. The initial 
magnetic field, which lies within the computational plane, is perpendicular to the cloud motion 
in the left panels (transverse case) and parallel to it in the right panels (aligned case). For each 
quantity and in each case two different times; namely t/tb c = 2, 6, are shown. 

Fig. 4b. — Same as in Fig. 4a but now for the oblique field case. The initial magnetic field makes 
an angle 6 = 45° with respect to the cloud motion. The resemblance between this case with the 
transverse case is noteworthy. 

Fig. 5a. — A light MHD cylindrical jet. The calculation has been done on a 256 x 1700 cylindrical 
grid with a computational domain r = [0, 1] and z = [0,6.64]. The sound speed of the ambient 
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medium, a am bient = 1, and its magnetic field is poloidal with (3 a mbient = 100. The jet has a radius 
of 30 cells, density contrast p jet I P ambient = 0.1, and Mach number Mj et = 20. The jet magnetic 
field is helical with maximum 0j et = 20 at the surface. The gray scale images show logarithmic gas 
density (upper frames) and logarithmic total magnetic pressure (lower frames) at t = 0.3, 0.8, 1.3 
1.8, and 2.2. White represents high values and black represents low values. 

Fig. 5b. — The same jet as in Fig. 5a. The arrows show velocity (upper frames) and r and z- 
components of magnetic field (lower frames) at t = 0.3, 0.8, 1.3 1.8, and 2.2. The length of velocity 
arrows is scaled as y/\v\ and that of magnetic field arrows as l-Bj 1 / 4 , in order to clarify the structures 
with small velocity and magnetic field magnitudes. 
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